Intrinsic localized modes in parametrically-driven arrays of nonlinear resonators 
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We study intrinsic localized modes (ILMs), or solitons, in arrays of parametrically-driven nonlinear 
resonators with application to microelectromechanical and nanoelectromechanical systems (MEMS 
and NEMS). The analysis is performed using an amplitude equation in the form of a nonlinear 
Schrodinger equation with a term corresponding to nonlinear damping (also known as a forced 
complex Ginzburg-Landau equation), which is derived directly from the underlying equations of 
motion of the coupled resonators, using the method of multiple scales. We investigate the creation, 
stability, and interaction of ILMs, show that they can form bound states, and that under certain 
conditions one ILM can split into two. Our findings are confirmed by simulations of the underlying 
equations of motion of the resonators, suggesting possible experimental tests of the theory. 

PACS numbers: 63.20.Pw, 05.45.-a, 62.25.-g, 85.85.+j 



I. INTRODUCTION 

The study of collective nonlinear dynamics of cou- 
pled mechanical resonators has been regaining atten- 
tion in recent years [l| thanks to advances in fabri- 
cation, transduction, and detection of microelectrome- 
chanical and nanoelectromechanical systems (MEMS and 
NEMS). Nonlinearity is readily observed in these sys- 
tems 0, i, i, i. El, 0, i, i, [13, [m, E, [il, and is even 
proposed as a way to detect quantum behavior [3, llH . 
Typical MEMS and NEMS resonators are characterized 
by extremely high frequencies — from hundreds of kHz to 
a few GHz [16|,[l3 — and relatively weak dissipation, with 
quality factors Q in the range of 10^ — 10^. For such de- 
vices, under external driving conditions, transients die 
out rapidly, making it is easy to acquire sufficient data 
to characterize the steady-state well. Because the basic 
physics of the individual elements is simple, and rele- 
vant parameters can readily be measured or calculated, 
the equations of motion describing the system can be es- 
tablished with confidence. This, and the fact that weak 
dissipation can be treated perturbatively, are a great ad- 
vantage for comparison between theory and experiment. 

Current technology enables the fabrication of large ar- 
rays, composed of hundreds or thousands of MEMS and 
NEMS devices, coupled by electric, magnetic, or elastic 
forces. These arrays offer new possibilities for quanti- 
tative studies of nonlinear dynamics in systems with an 
intermediate number of degrees of freedom — much larger 
than one can deal with in macroscopic experiments, yet 
much smaller than one confronts when considering non- 
linear aspects of phonon dynamics in a crystal. Our stud- 
ies of collective nonlinear dynamics of MEMS and NEMS 
were originally motivated by the experiment of Buks and 
Roukes [18|, in which an array of 67 doubly-clamped mi- 
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cromechanical gold beams was parametrically excited by 
modulating the strength of an externally-controlled elec- 
trostatic coupling between neighboring beams. These 
studies have led to a quantitative understanding of the 
collective response of such an array, providing explicit bi- 
furcation diagrams that explain the transitions between 
different extended modes of the array as the strength 
and frequency of the external drive are varied quasistat- 
ically [1^, [20] . We have also considered more general is- 
sues such as the nonlinear competition between extended 
modes, or patterns, of the system — when many such pat- 
terns are simultaneously stable — as the external driving 
parameters are changed abruptly or ramped as a func- 
tion of time Furthermore, we have investigated the 
synchronization that may occur in coupled arrays of non- 
identical nonlinear resonators, based on the ability of 
nonlinear resonators to tune their frequency by chang- 
ing their oscillation amplitude 0, . 

Here we focus on a different type of nonlinear states, 
namely, intrinsic localized modes (ILMs), also known as 
discrete breathers or lattice solitons [3 [25|, [2^ . These 
localized states are intrinsic in the sense that they arise 
from the inherent nonlinearity of the resonators, rather 
than from extrinsically-imposed disorder as in the case 
of Anderson localization. ILMs have been observed by 
Sato et al [13, [H, [H, [sO, [111, [32[ in driven arrays of 
micromechanical resonators. They have also been ob- 
served in a wide range of other physical systems including 
coupled arrays of Josephson junctions [H, [sj , coupled 
optical waveguides [35I. [36l . W\ , two-dimensional nonlin- 
ear photonic crystals [38|, higlily- nonlinear atomic lat- 
tices [39I, and antiferromagnets [40,|4l|. Thus, the ability 
to perform a quantitative comparison between our theory 
and future experiments with large arrays of MEMS and 
NEMS resonators, may have consequences far beyond the 
framework of mechanical systems considered here. 

We aim to predict the actual physical parameters, 
in realistic arrays of MEMS and NEMS resonators, for 
which ILMs can form and sustain themselves. Such pre- 
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dictions may have practical consequences for actual ap- 
plications exploiting self-localization to focus energy, and 
others that may want to avoid energy focusing, for exam- 
ple in cases where very large oscillation amplitudes may 
lead to mechanical failure. Although quantitative anal- 
ysis can be carried out directly in the framework of the 
underlying oscillator equations of motion, it is instructive 
to formulate the analysis in terms of an amplitude equa- 
tion, as done previously for extended modes [l|, IIO, |2l|. 
This allows one to display the range of stable ILMs on 
a reduced diagram, helping to describe the general qual- 
itative behavior as physical parameters are varied. In 
Sec. ini we describe the derivation of such a single am- 
plitude equation from the coupled equations of motion 
that model an array of nonlinear resonators. This ampli- 
tude equation is obtained in the form of a parametrically- 
driven damped nonlinear Schrodinger equation with an 
additional nonlinear damping term, also known as the 
forced complex Ginzburg-Landau equation. In most 
physical systems, the dissipation of energy is modeled by 
a linear damping term. However, it has been established, 
both theoretically [l|, [l^ and experimentally [42[, that 
nonlinear damping is important for correct modeling of 
certain high-Q nonlinear MEMS and NEMS resonators. 

In Sec. mil we argue that exact soliton solutions that 
exist in the absence of nonlinear damping can be contin- 
ued to solve the full amplitude equation, with nonlinear 
damping. We derive an approximate analytical expres- 
sion for these solitons, and then use it to find the exact 
soliton solutions numerically. A linear stability analysis 
of these soliton solutions follows in Sec.[lVl In Sec.lVlwe 
consider the dynamical formation of solitons and study 
the effects of nonlinear damping on the modulational in- 
stability of non-zero uniform solutions of the amplitude 
equation. In Sec. [Vl] we study the interactions between 
pairs of solitons; and in Sec. I VIII demonstrate the split- 
ting of a single soliton into two separate ones. Finally, 
we show in Sec. IVIIII that solitons of the full amplitude 
equation can form stable bound states. As emphasized 
in the concluding remarks in Sec. \IX\ all the phenom- 
ena demonstrated through the analysis of the amplitude 
equation are accurately reproduced in simulations of the 
underlying equations of motion of the coupled resonators, 
and therefore should also be reproducible in actual exper- 
iments with MEMS and NEMS arrays. 



II. DERIVATION OF THE AMPLITUDE 
EQUATION 

Lifshitz and Cross [l^ modeled the array of cou- 
pled nonlinear resonators that was studied by Buks and 



Roukes (iS^ using the equations of motion 

- {Un - Un-lf{Un - Un-l)] = 0, (1) 

where Un describes the deviation of the n^^ resonator 
from its equilibrium position, with n = 1 . . . N ^ and fixed 
boundary conditions uq = un-\-i = 0. Detailed argu- 
ments for the particular choice of terms introduced into 
the equations of motion are discussed in Ref. [l^. The 
terms include an elastic restoring force with both lin- 
ear and cubic contributions (whose coefficients are both 
scaled to 1), a dc electrostatic nearest-neighbor cou- 
pling term with a small ac component responsible for 
the parametric excitation (with coefficients D and H 
respectively), and linear as well as cubic nonlinear dis- 
sipation terms. The nonlinear elastic term is positive, 
indicating a stiffening of the resonators with increasing 
displacement, which is the common situation when using 
doubly-clamped beams. Both dissipation terms are taken 
in the nearest-neighbor form, which is motivated by the 
experimental indication that most of the dissipation orig- 
inates from the electrostatic interaction between adjacent 
beams. Note that the electrostatic attractive force, act- 
ing between neighboring beams, decays with distance, 
and thus acts to soften the elastic restoring force. For 
this reason the sign in front of the coupling coefficient 
D is positive, and accordingly the dispersion curve in 
the linear regime features a negative slope, or a negative 
group velocity. 

In more recent implementations the electric cur- 
rent damping has been reduced, and the parametric drive 
is applied piezoelectrically directly to each resonator, 
simplifying the equations modeling the array, 

Un + Q~^Un + [l - H CO^{2uJpt)\Un + + WnUn 

- ^D{Un^l-2Un^Un-l) =0. (2) 

The negative sign before the coupling coefficient D mod- 
els elastic coupling between adjacent beams, which is 
stronger as the separation between neighbors increases, 
thus acting to stiffen the resonators. Here the dis- 
persion curve has a positive slope, or a positive group 
velocity. The coupling mechanism in the experimen- 
tal setups in which ILMs have been observed is of this 

kind [23, m, m, [33, [m, [31 . 

Because the quality factor Q of a typical MEMS or 
NEMS resonator is high we follow the practice [l|, [H, 
of using it to define a small expansion parameter = 
67, with e <C 1, and 7 of order unity. The driving ampli- 
tude is then expressed as = e/i, with h of order unity, 
in anticipation of the fact that parametric oscillations 
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at half the driving frequency require a driving amphtude 
which is of the same order as the Unear damping rate |43|] . 

An experimental protocol for producing ILMs in an 
array of resonators with a stiffening nonlinearity — albeit 
not the one we use below — is to drive the array at the 
highest-frequency extended mode. As the resonators are 
collectively oscillating at this mode, the frequency is 
raised further which results in an increase of the oscil- 
lation amplitude up to a point in which the extended 
pattern breaks into localized modes [23, [lO] • With this 
in mind — and concentrating on the case of elastic cou- 
pling where the highest-frequency mode uo = Vl + 21) is 
the staggered mode, in which adjacent resonators oscil- 
late out of phase — we write the displacement of the n^^ 
resonator as 



with slow temporal and spatial variables T = et and 
Xn — ^^^Ti^ and c.c. standing for the complex conju- 
gate expression. We take the parametric drive frequency 
to be close to twice (jJ by setting cc;^ = cc; + el^/2, intro- 
duce a continuous spatial variable X in place of X^, and 
substitute the ansatz ([3]) into the equations of motion ([2]) 
term by term. Up to order e"!^ we have 



.1/2 



dT) 



(4a) 



.1/2 



c.c. 



1/2^ + £ ^i^^t-i^n) 

dX 2 I 



(4b) 



c.c, 



e7u„ = e^/^^iw^e*''^*-™) + c.c, 
^3 =e3/23|^/,|2^e^M-™)+0(( 



(4c) 
(4d) 

)+cc, (4e) 
(4f) 



where 0(6 ^^ , e ) are fast oscillating terms with tem- 
poral frequency 3uj or spatial wavenumber Sir. 

At order e^/^ the equations of motion ([2|) are satisfied 
trivially. However, at order e^/^^ q^iq must apply a solv- 
ability condition [44] , requiring all terms proportional to 
gz(a;^-7rn) vanish. It is this condition that leads to a 
partial differential equation (PDE) describing the slow 
dynamics of the amplitudes of the resonators, 

(5) 



Note that while e^^^*+^^) = e^^^^-^^), if we were to con- 
sider an arbitrary mode of wave number k instead of tt, 
the parametric term would have forced us to apply an- 
other solvability condition, requiring terms proportional 
to 6*^^^^+^^) to vanish. In that case an ansatz based on 
counter propagating waves is considered as the 0{e^^'^) 
solution for Un and a system of two coupled amplitude 
equations emerges, as shown in Ref. |2Qi] . 
By means of rescaling. 
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D ^ 2 
Y T = T 

2un ' n ' 



h = 2u;Qh^ 7 = 1^7, fj 
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we transform Eq. ([5]) into a normalized form. 



(3) 

dT 



■ i-yip - (2 + ir])\ipfip + hip*e 



2iT 



(6) 



(7) 



We then perform one final transformation ipc'^^ and 
arrive at an autonomous PDE, which is the amplitude 
equation that we study in the remainder of this work. 



dT 



dX^ 



-^{l-i-f)^-{2^ir])\^\^^^h^\ (8) 



Equation ([7j) with 77 = is called the parametrically- 
driven damped nonlinear Schrodinger equation 
(PDNLS). It models parametrically driven media 
in hydrodynamics [45^, 46, 47_, ^8] and optics [H, [s^ , and 
was also used as an amplitude equation to study localized 
structures in arrays of coupled pendulums [5l|, (HH, [s^. 
Recently, a pair of linearly-coupled PDNLS equations 
was used to model coupled dual-core wave guides [H^ . 
Eq. (|8]) has the form of a forced complex Ginzburg- 
Landau equation [55[ but with specific coefficients that 
are derived, via the scaling performed in (|3]) and (|6]), 
from the underlying equations of motion ([2]). 

We note that considering the equations of motion ([1]) 
(yet still with a negative sign before D) instead of Eqs. ([2]) 
leads to the same equation (j5j) as above, but with differ- 
ent coefficients (a factor of 2 multiplying h and 7, and a 
factor of 8 multiplying fj). Thus, applying modified scal- 
ing ([6]) yields exactly the same amplitude equation ([8|). 



III. SOLITONS IN THE PRESENCE OF 
NONLINEAR DAMPING 

A. Continuation of the PDNLS solitons to 77 > 

A remarkable feature the amplitude equation ([8|) is 
that for ?7 = it has exact time-independent solitonic 
solutions, as shown by Barashenkov et al. [561 ], 

^±{X) = A±e-^®=^sech [A± {X - Xq)] , 



(9) 



where Xq is an arbitrary position of the soliton, and 



Al = i± ^/W^, cos(2e±) = ±^/i - 



(10) 
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This pair of solitonic solutions exists for 7 < /i. It was 
shown in [s^ that the ^_ sohton is unstable for all values 
of 7 and while the ^+ soliton is stable in a certain 
parameter range. A simple linear stability analysis shows 
that the zero solution ip{X) = 0, which exists for all 
parameter values, is stable only for h < y^l + 7^. This 
inequality also determines an upper stability limit for 
localized solutions of Eq. ([8]) that decay exponentially to 
zero on either side. 

The aim of this section is to show that the PDNLS soli- 
tons ^± can be continued to nonzero nonlinear damping 
77. A similar calculation was performed by Barashenkov 
et al. where they considered the addition of a spectral 
filtering term —icd'^il;/dX'^ to the PDNLS. We do so by 
expanding the stationary solutions of the full amplitude 
equation ([8|) in powers of 7^, which is assumed small, to 
get 



(11) 



SO that (j)o = |^±|- Denoting (j)n = Un -\- ivn^ with real 
Un and Vn, substituting i/j into Eq. ([8]), and comparing 
powers of r] yields equations of the form 



where 



L± 



_f-d%-6ul + Al 
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-d% - 2ul + A 



(12) 



(13) 



One can use Eq. (p!2|) iteratively to find the n^^-order 
correction given all lower-order ones, provided that 
the right-hand side is orthogonal to the null subspace, 
or kernel, of the adjoint operator L^, if such a sub- 
space exists. Indeed, in the relevant parameter range, 
7 < < + 7^, the adjoint operator has one zero 
eigenvalue, and the corresponding eigenvector consists of 
only odd functions of X — Xq [13, [H^. On the other hand, 
one can verify that the functions Fn and Gn , which orig- 
inate from the nonlinear terms in Eq. ([8]), include only 
natural powers of iio, '^0, • • • , Un-i^Vn-i- Therefore, and 
because uq = |^±|, = 0, and L± is parity preserv- 
ing, the right hand side of Eq. ([12]) consists of only even 
functions of X — Xq, for any n. This suggests that it 
is possible to continue the PDNLS solitons ^± to solve 
Eq. ([8]) up to any order in r]. This can also be done in 
practice, by calculating Fn and Gn symbolically, express- 
ing L± as a discrete matrix, and inverting it to find Un 
and Vn in each iteration, although we do not follow this 
procedure here. 



B. Approximate analytical solitons with 77 > 

Motivated by the arguments above, we wish to con- 
struct an approximate analytical expression for the lo- 
calized solution of the full amplitude equation JHl) , im- 
plementing the method of Barashenkov et al. [53] • To 



this end, we consider a function of the same form as 

V^(X, T) = a(T)e-^^(^)sech [a(T) (X - Xq)] , (14) 

except that a and 6 are now time-dependent. We multi- 
ply Eq. ([8]) by ?/^*, subtract the complex conjugate of the 
resulting equation and get 



]+h[{rf-i^^] 



(15) 



By substituting i/j = l^/^le"*^, integrating over X' = X — 
Xq, and assuming that i/j ^ and dilj/dX ^ as |X| 
00, we obtain a spatially-independent integral equation 



d 

dT 



\ilj\^dX' = 2 J \ij\^[hsm{2x)--f]dX'-2r] J |V^|'^dX^ 

(16) 

Substituting the ansatz ([T4|) into Eq. (fT6|) . we obtain the 
time evolution equation for a 



da 
dT 



(17) 



where fj = 2r]/3. The time evolution equation for 6 is 
derived in a similar way by multiplying Eq. ([8|) by ?/^*, 
adding the complex conjugate of the resulting equation, 
substituting the ansatz ([T4|) . and integrating over space 
to yield 



df) 

— = hcos{29) + l 



(18) 



Equations ([T7|) and (p^ have the same form as the 
equations obtained in [57|, whose fixed points are 



1 - 7^ ± V/l2(l+7?2) - (7 + ^)^ 



1 



which has to be positive, and 



/icos(26>±) 
h^\n{20±) 



fj2 



(19) 



1, 



7 + ^tt±. 



(20) 



A linear analysis of these stationary points shows that 
(a+,6>+) and (a ,0-) are a stable node and a saddle, 
respectively [S^- The saddle-node bifurcation point of 
these solutions occurs at 



hsn{ri) 



7 + ^ 



- fp' 



where fj = -77, 
o 



(21) 



as long as 77) < 1. This is the approximate minimal driv- 
ing strength required to support a localized structure in 
the array, in the presence of linear and nonlinear dissipa- 
tion. 

The approximate stable localized solution of the am- 
plitude equation ([8]) is therefore given by 

^app(X) = a+e-^^+sech(a+(X - Xo)). (22) 
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Substituting this expression into Eq. (|3]) yields an ap- 
proximate expression for the displacements of the actual 
resonators in the array, 



Un{t) 



2\l — - — a+sech 



X cos {cOpt — irn 




n- Xo 



(23) 



C. Numerical solutions for solitons with 77 > 

To obtain accurate solutions we solve the amplitude 
equation, as well as the underlying discrete equations of 
motion, numerically. The equations of motion ([2]) are ini- 
tiated with the approximate expression ([23|) at a value of 
h just above the saddle node hsn ([2T]) . We then perform 
a quasistatic upward sweep of /i, raising h in small incre- 
ments and waiting for transients to decay at each step. 
To obtain the stationary solution of the amplitude equa- 
tion we set dip/dT = in Eq. ([8]) and solve it numerically 
as a boundary value problem over an interval of length L, 
with boundary conditions ?/^(X = 0) = ^^(X = L) = 0. 
We use the approximate expression ?/^app(-^) [Eq. 
as an initial guess. 

In Fig. [1] we compare the integral measure 



Hfix)} 



2 7-, 



f\X)dX 



(24) 



for the different localized solutions, where for the analyti- 
cal solutions f{X) = asech(aX) cos(^) and / = acos^(6>), 
whereas for the stationary numerical solution V^(X) of 
the amplitude equation ([8]) /(X) = Re?/;(X), and for the 
numerical steady-state solution of the equations of mo- 
tion ([2]) /(X) is the appropriately scaled magnitude \un\ 
of the n^^ resonator measured at times tm = ^nm/up. 
Fig. [1] shows good agreement between the numerical solu- 
tion of the equations of motion ([2]) , and the approximate 
and numerical solutions of the amplitude equation (|8]). 

Fig. [2] shows the phase of the numerical solution of the 
amplitude equation (|8]), calculated as 



arctan 



Im?/^\ 



(25) 



One can see that while the approximate phase 6>+ pro- 
vides a good estimate for the actual phase of the solu- 
tion near the peak of the soliton, the phase asymptotes 
to 6- = 7r/2 — 6+ as the soliton's amplitude decays to 
zero. This is a surprising result since it might have been 
expected that the phase of a numerical continuation of 
the ^+ solution would tend back to 9+ as the amplitude 
drops to zero, eliminating the nonlinear damping. How- 
ever, a similar situation was observed in a bound state of 
two ^+ solitons in the PDNLS equation [5^. 




FIG. 1: (Color online) The integral measure /{/(X)} as a 
function of h. Red outer solid and dashed lines represent the 
exact analytical solutions of the PDNLS equation without 
nonlinear damping cos^ G+ and A- cos^ G-, respectively. 
Black inner solid and dashed lines represent the approximate 
analytical solutions with nonlinear damping a+ cos^ 0^ and 
a-cos^^_, respectively. These analytical lines end at the 
upper stability boundary h = \/l -\- 7^. The points desig- 
nated by crosses and circles (+ and o) are taken, respectively, 
from the stationary numerical solution of the amplitude equa- 
tion JS]), and from the numerical solution of the equations 
of motion ([2]), as elaborated in the text. The inset shows 
the absolute value of the profile of the solution for h = 0.87 
where os are results of the numerical solution of the equa- 
tions of motion (|2)). The solid line shows the real part of the 
numerical solution of the amplitude equation dS]), scaled 
by a factor of 2y^2eLjQ/3. The scaled analytical approxi- 
mation (|22p is indistinguishable from the solid line in this 
plot. The dot dashed line shows the scaled analytical solu- 
tion (|9|) in the absence of nonlinear damping. The parameters 
are 7 0.5, D 0.25, e 0.01, cjp 1.002cj,ry 0.1, and 
N = 399. 



IV. LINEAR STABILITY ANALYSIS OF 
SOLITONS 



Having identified an upper stability boundary h = 
Y^l + 7^ and an approximate lower existence boundary, 
given by Eq. ([2T]) . we turn to examine the stability of 
the localized solution within these boundaries. For this 
purpose, we substitute into Eq. ([8|) ^^(X, T) = V^(X) + 
SiIj{X^ T), where V^(X) could be any steady-state solution 
of the equation — in this case the stationary localized so- 
lution, which is obtained numerically — and ^'^(X, T) is a 
small perturbation. We linearize in ^?/^(X, T), substitut- 
ing i/j = R -\- il and Si/j = U -\- iV with real /, [/, and 
and obtain the equation 



dT \V 



H 



(26) 
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FIG. 2: (Color online) The phase of the numerical solution 
of the amplitude equation (|8]), and the analytical exact and 
approximate phases, as indicated in the legend. Near the peak 
of the soliton, the phase of the numerical solution is close to 
^+ and it asymptotes to G- as the soliton's amplitude decays 
to zero. Parameters are the same as in the inset of Fig. [1] 
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FIG. 3: (Color online) Stability diagram for localized solu- 
tions of the amplitude equation Q in the /i vs. 7 plane. The 
dotted line is the lower existence boundary for ?7 = 0, namely 
h — ^. The dashed-dotted line is the approximate low bound- 
ary for Tj = 0.1, given by Eq. (|2T]) . Above the solid line the 
^+ solution of the PDNLS equation with 77 = is unstable 
with respect to a Hopf bifurcation [56|] . The dashed line is the 
line h = -^/l + 7^ above which the zero solution is unstable. 
Red *s are points for which the matrix J~^H has a pair of 
complex conjugate eigenvalues with a positive real part for 
T] = 0.1, hence the soliton solution ip{X) is unstable. Blue 
dots represent points for which the solution '0(X) is stable 
according to the linear analysis. The four black circles la- 
beled (a)-(d) indicate parameter values corresponding to the 
numerical simulations of the equations of motion (|2]) shown 
in Figs.|4i;a)-|4i;d). 
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FIG. 4: (Color online) Results of the numerical solution of the 
equations of motion (|2]) for 77 = 0.1, 7 = 0.1225, and values 
of h labeled as (a)-(d) in Fig. O The solutions are plotted at 
times tm = 27rm/a;p, with integer values m shown on the ver- 
tical axis, (a) h = 0.17 is below the approximate low bound- 
ary ([2T]) but above the low boundary for 77 = 0. One sees 
that the localized structure decays to zero, (b) h = 0.1988 is 
above the approximate low boundary, where linear stability 
analysis predicts that the soliton is stable, (c) For h = 0.35 
the stationary soliton is unstable, and an oscillating localized 
solution is formed instead. In (d) h = 0.85 and the soliton 
is stable again. The stability is due to nonlinear damping, 
without which the soliton is unstable as demonstrated in (e) 
where h — 0.85 and ry = 0. All unspecified parameters are the 
same as in Fig. [1] 



where 

Hii = -dx - 6i?2 - 2/2 + 1 + /i + 2riRI, 

Hu = + 3/2) - 4i?/ + 7, 

H21 = -r?(3/?2 + /2) - 4/?/ - 7, 

H22 = -d\-QI^ - 2R^ + l-h-2r]RI. (28) 

By expressing the small perturbations as U {X, T) = 
Re[u{X)e^^] and V{X,T) = Re[v{X)e^^], where \ u, 
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and V are complex, we arrive at the eigenvalue problem 




When is obtained numerically, the eigenvalues of 

the matrix J~^H, describing the growth of perturba- 
tions T), are found by performing a spatial dis- 
cretization of Eq. ([29|) and diagonalizing J~^H numeri- 
cally (see [60| for details). The stability diagram of both 
the analytical solution ^+ for 77 = [56| and the numer- 
ical solution for r] = 0.1 are displayed in Fig. [3l 
These results are verified at a few points by numerical 
integration of the equations of motion ([2]), as shown in 
Fig. [2 

Fig. [3] highlights the effects of nonlinear damping on 
localized solutions. The first effect is to raise the lower 
existence boundary. This is explained by the fact that the 
additional energy lost through nonlinear damping has to 
be compensated by an increase in the strength of the 
parametric drive, as predicted by the approximate ex- 
pression ([2T]) . The second effect is that nonlinear damp- 
ing increases the area in the (/i, 7) parameter space where 
solitons are stable (blue dots). In particular, the shape 
of the unstable region for 77 > (red *s) becomes qual- 
itatively different. There are values of 7 for which an 
increase in the drive amplitude h initially induces an in- 
stability of the soliton, while upon further increase of h 
the soliton regains its stability. This can be explained 
by noting that the amplitude of the soliton — given ap- 
proximately by Eq. (fT9|) — increases as h becomes larger, 
thereby enhancing the effect of nonlinear damping. This 
increase of damping exerts a similar stabilizing effect as 
that of increasing 7 in the absence of nonlinear damping. 

The effect of regained stability with the increase of h 
does not occur in the absence of nonlinear damping, as 
the ^+ soliton is unstable for all parameter values above 
the solid black line in Fig. [3] [56|. Different solutions 
of the PDNLS equation above this instability threshold 
were found by Bondila et al. [6lj to be localized solutions 
that oscillate in time with different periods and chaotic 
solutions, in addition to the zero solution. As we in- 
crease the nonlinear damping coefficient 77, the region in 
(/i, 7) space, in which the single soliton becomes unstable 
against these alternative solutions, shrinks in size. 



V. DYNAMICAL FORMATION OF SOLITONS 

A. Self-trapping of solitons 

It is not obvious how dynamically to form solitons 
starting with a motionless array of resonators, as one 
needs to take the system sufficiently far from the basin 
of attraction of the zero solution V^(X) = 0, which is also 
stable whenever solitons are stable. The most direct pro- 
cedure for avoiding the zero solution, starting from weak 
random noise, is to drive the system with h > yl + 7^. 




50 100 150 200 50 100 150 200 



n n 

FIG. 5: (Color online) Numerical simulation of the coupled 
equations of motion (|2]) showing the dynamical creation of 
solitons. Linear damping is set to 7 = 1 and nonlinear damp- 
ing to ?7 = 0.3. All other parameters are the same as in Fig.[T] 
Plotted are the absolute values of the displacements of the res- 
onators, which alternate between positive and negative values. 
Left panels show the complete time evolution, with m count- 
ing the number of drive periods. Right panels show the initial 
(black dots) and final (blue circles) states along with the ana- 
lytical form of the solitons (green solid line), using only their 
central positions Xq as fitting parameters. Top panels: A 
simulation of 199 resonators with fixed boundary conditions 
is initiated with random noise and a drive ampli tude of /i = 5, 
which is above the upper stability limit, h = ^^1+7^ = ^/2, 
for both the zero-state and the solitons. At time m — 600 
drive periods, after some non-zero transient (black +s in the 
right panel) has developed, the drive amplitude is lowered 
to h = 1.35 < v^, yielding stable solitons. Bottom panels: 
A simulation of 200 resonators with periodic boundary con- 
ditions is initiated with the uniform non-zero solution and 
a drive amplitude of h = 1.3, which is above the stability 
threshold §6^, hth ^ 1.26, for this state. After m = 10000 
drive periods during which the uniform state remains stable, 
the drive amplitude is lowered to /i = 1.2 < hth, yielding 
stable solitons. 



so neither the zero solution nor the soliton solutions are 
stable. As a consequence, a non-zero pattern develops. 
Stable solitons can then be formed by lowering the drive 
amplitude to a value h < y^l + 7^ for which the zero 
solution and the soliton solutions are both stable, if the 
non-zero pattern that was obtained is outside the basin 
of attraction of the zero solution. 

This simple procedure — which could be implemented 
experimentally in a straightforward manner — is demon- 
strated in the top panels of Fig. \5\ showing a numeri- 
cal simulation of the equations of motion ([2]) with fixed 
boundary conditions, using N = 199 resonators. One 
can see that the initial transient that forms becomes un- 
stable upon lowering the drive amplitude, giving rise to 
the formation of a number of solitons. Note that be- 
fore reaching steady state a pair of solitons merges into 
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one, and another pair attracts and forms a bound state. 
Both of these effects are studied below. The emerging 
isolated solitons agree well with the approximate analyt- 
ical form ([23|) . determined earlier, with only their central 
positions Xq used as fitting parameters. 



B. Modulational instability of uniform states 

A more controlled procedure for generating solitons 
would be to initiate the array in a particular non-zero 
state and drive it outside its known stability boundaries. 
This has been considered in the past in systems without 
nonlinear dam ping, us in^ the non-zero uniform solution 
of the PDNLS [571 1621. [63l. [63|. However, it is known for 
systems with 77 = that the uniform solution is always 
unstable against weak modulations and so may be dif- 
ficult to access dynamically. We wish to examine here 
whether the non-zero uniform solution may be stabilized 
with the help of nonlinear damping (77 7^ 0), thereby mak- 
ing it accessible dynamically and possibly opening an ad- 
ditional experimental route to the formation of solitons. 

Indeed, the amplitude equation (|8|) admits a pair of 
non-zero spatially-uniform solutions of the form 



ip± = a±e 



(30) 



Substituted into the perturbative expansion (|3]), this 
yields an oscillation of the array in its staggered mode 
with wave number tt, about which we initially expanded 
our solution. If we impose fixed boundary conditions the 
staggered mode will be modified near the boundaries to 
accommodate these conditions, but would otherwise re- 
main unchanged in the bulk of the system. 

Letting fj = r]/2 and substituting the uniform solution 
(|3Q|) into the amplitude equation (|8]) yields 



2a^ 



1 - 7^ ± ^/h^O^TWh^WTW 



1 + 7^^ 

which has to be positive, and 



(31) 



hcos{2e±) 
hsm{2e±) 



2a\ 
7 + 



-1, 
r?(24). 



(32) 



For 77^ < 1 both solutions exist, and a saddle- 
node bifurcation — obtained by setting the square-root in 
Eq. ([3T]) to zero — occurs at 



hsniv) 



7 + ^ 



77^ 



where ^ ~ 2* ^^^^ 



For 77) > 1 the bifurcation from the zero solution be- 
comes supercritical, occurring at the instability boundary 
of the zero solution h = y^l + 7^. Note that apart from 
rescaling 7^ by a factor of 3/4 and a± by a factor of >/2 
these expressions are identical to those for the approxi- 
mate amplitude and phase of the soliton solutions ([19]- 
[2ID. 




FIG. 6: Stability boundary of the large- amplitude uniform 
solution -0+ for 7 = 1. The solution exists above the 
dashed and dot-dashed curves. The dashed curve is the 
saddle-node hsn{ff) [Eq. (|33)) ]. which is replaced at the point 
marked with a small square by the dot-dashed curve, indi- 
cating the supercritical bifurcation from the zero solution at 
h = + = ./2. The solid curve shows hth{v) [Eq- (|36p ] 
above which V < 0. It coincides with the dashed curve at 
fjc — —7 + + 72 — — 1, indicated by a small circle. 
In the dark-gray region -0+ is stable because both zeros of 
the quadratic function in (|34p are complex. In the light-gray 
region -0+ is stable because both zeros are real and negative. 
This implies that for fj > fjc the large- amplitude solution is 
always stable, and for fj < fjc the drive h has to exceed the 
threshold value hth{fj) [Eq. (|36p ] before the solution becomes 
stable. Recall that fj — rj/2. 



The modulational instability of the uniform solutions 
can be evaluated [s^ by adding perturbations of the form 
exp(±iA:X) and calculating their growth rates using the 
eigenvalues of the matrix J~^H^ obtained from Eqs. ([27|) 
and ([28|) by substituting —d\ = k^^R = a± cos^±, and 
I = —a± sin(9±. The uniform solutions are stable against 
such modulations as long as the larger of the real parts of 
the two eigenvalues of J~^H is not positive for any real k. 
This requirement translates to satisfying the inequality 



s^ + 2(l-4a^)5±8a^V^2(l + 7^2) _ ^ > q (34) 

for any non- negative s = k'^, where a positive sign is 
assumed for the square root which is real for h > hsniv)- 

As expected, the small-amplitude uniform solution 
is unstable even against a uniform perturbation, as the 
left-hand side of (|34|) is negative for /c = 0. For the 
large-amplitude uniform solution ^/S^ the inequality (|34|) 
is satisfied if either both zeros of the quadratic function 
of s on the left-hand side are complex, or both are real 
and non-positive. The first condition is satisfied if the 
discriminant 



V = 4- 32f]al (7 + 2f]al) 



(35) 



is negative, and the second condition is satisfied if P > 
and 4a^ < 1 (because the constant term in the quadratic 
function is positive). Clearly, for 7^ = these stability 
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conditions are not satisfied, and the large- amplitude uni- 
form solution is modulationally unstable, in agreement 
with known results [53,[62[. However, we indeed find that 
?/^+ can be stabilized with the help of nonlinear damping, 
as shown in Fig. [6l If fj > fjc = —7 + y^l + 7^ then tjj-^ 
is stable everywhere. For weaker nonlinear damping the 
drive h must exceed a threshold value 



^|l + 272 (1 + ^^) +477^ + 57^^ 

. 1/2 

-2 (7 + 27/ - 7f ) + , 



(36) 



determined by substituting the expression for 
[Eq. (|3T]) ] into the discriminant (|35|) . and setting V = 0. 

If the inequality ([34j) is not satisfied, the uniform state 
is modulationally unstable, and the modulation whose 
growth rate is fastest is expected to appear. This mod- 
ulation corresponds to the minimum of the quadratic 
function on the left-hand side of (|34|) . with wave num- 
ber kfast = \J^o\ - 1- 

We demonstrate the use of the stable uniform solution 
in the dynamical formation of solitons in the bottom pan- 
els of Fig. [5l showing a numerical simulation of the equa- 
tions of motion ([2]) with periodic boundary conditions, 
using N = 200 resonators. The array is initiated with the 
large-amplitude uniform solution and is driven within the 
stability boundary of this state. After a long time during 
which the uniform solution remains stable, the drive am- 
plitude is lowered below the stability threshold (|36|) for 
this solution, but within the stability boundaries of the 
soliton solutions, and solitons are formed via a modula- 
tion of the unstable uniform state. The wavenumber of 
the modulation that is observed numerically agrees with 
the predicted value kfast to within rounding to the near- 
est mode satisfying the periodic boundary conditions. 



VI. SOLITON INTERACTIONS 

After finding a family of stable soliton solutions of 
Eq. ([8]), it is natural to consider the interaction between 
them. Soliton interactions were studied in detail in the 
integrable NLS equation, corresponding to7 = /i = 7/ = 
in Eq. ([7j) [H, [6^. It was found that the interaction 
of initially stationary solitons depends on their relative 
phase, with in-phase and out-of-phase solitons attract- 
ing and repelling each other, respectively. This property 
is generic and is not predicated on the integrability of 
the underlying equations. It is also valid for multidimen- 
sional equations [67]. In the presence of additional effects 
such as amplification and damping the soliton interac- 
tion problem is not amenable to a complete analytical 
study [H, [1^. However, it is possible to analyze the in- 
teraction between two weakly overlapping ^+ solitons by 
regarding the overlapping nonlinear terms — arising from 
the substitution of a two-soliton solution into the PDNLS 




(a) 





(c) 



FIG. 7: (Color online) Numerical simulations of soliton in- 
teraction. Left and right panels display results obtained, re- 
spectively, from numerical simulations of the amplitude equa- 
tion JS]) and of the underlying equations of motion (|2]), with 
h = 0.7 and the remaining parameters as in Fig. [1] (a) 
and (b): Attraction between two in-phase solitons and their 
merger into a single one, after half the energy has been dissi- 
pated, (c) and (d): Repulsion between out-of-phase solitons. 



equation — as small perturbations [49|, |70|, lDJ. In this 
case as well, in-phase solitons attract each other, whereas 
out-of-phase solitons repel. Using equations derived in 
Refs. [sil, m, [7l| it is easy to show that adding a small 
nonlinear damping term does not induce any motion on 
the solitons, hence one may expect to see the same type of 
phase-dependent interaction in the full amplitude equa- 
tion dH) with 7/ > 0. 

This is indeed verified, as shown in Fig. [71 which 
presents the results of a numerical integration of the 
equations of motion ([2|), and of the amplitude equa- 
tion dHl), simulated as a PDE with initial conditions 



^± = ^app ( X - (1 - r)^] ± V^app fx - (1 + r)^ 



(37) 

where L = {N -\-l){2ujne/Dy^'^ is the scaled length of an 
array of resonators, and rL is the distance between the 
centers of the solitons. Note that one time unit T = 1 of 
the amplitude equation is equal to u;p/{27T{cOp — co)) pe- 
riods of parametric oscillations in the original equations 
of motion. For the parameters used throughout this pa- 
per Ljp/{27T{cOp — co)) ^ 80, and one can verify that this 
is approximately the ratio between the vertical axes of 
the solutions of the amplitude equation [Fig. [TJa) and 

(c) ] and those of the equations of motion [Fig. [71(b) and 

(d) ]. Also note that the ratio between the peak heights of 
the soliton solutions of the amplitude equation and those 
of the equations of motion in Fig. [3 is approximately 
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FIG. 8: (Color online) Annihilation of a pair of strongly over- 
lapping solitons. Left and right panels display results ob- 
tained, respectively, from numerical simulations of the am- 
plitude equation dS]) and of the underlying equations of mo- 
tion ([2|), initiated with a separation given by the indicated 
values of rc above which annihilation was not observed, (a) 
and (b): An initial attraction followed by the annihilation of 
a pair of in-phase solitons. (c) and (d): An initial repulsion 
followed by the annihilation of a pair of out-of-phase solitons. 
Parameters are the same as in Fig[71 



^/e = 0.1, as expected from Eq. (|3]). 

Another effect, which is demonstrated in Fig. El is that 
if the solitons strongly overlap, and r is smaller than 
some critical distance Tc, they annihilate into the zero 
state, for both the in-phase and out-of-phase pairs. For 
the parameters of Fig. El Tc — 0.021. The annihilation 
of the out-of-phase pair is easily understood because the 
strongly overlapping solitons cancel each other. However, 
the mutual destruction of the in-phase pair is a less ob- 
vious effect, indicating that for some reason the initial 
conditions for r < Vc are in the basin of attraction of the 
zero solution, and not in that of the stable single soliton 
solution, contrary to the case oir > 



VII. SPLITTING SOLITONS 

The Galilean invariance of the NLS equation ad- 
mits the motion of any solution at a constant veloc- 
ity. The parametric drive hip'' breaks this property of 
the equation. Nevertheless, stable traveling solitons in 
the parametrically-driven (but undamped) NLS equation 
were obtained in a numerical form by Barashenkov et 
al. [72[. We have attempted to do the same with solu- 
tions of the full amplitude equation ([8|) by multiplying 
the approximate solution V^app by g-^^C^-^o)^ thereby 
boosting it. We have concluded that such a boost may 
set the soliton into transient motion, but eventually it 



FIG. 9: (Color online) Simulations initiated with a boosted 
soliton. Left and right panels display results obtained, re- 
spectively, from numerical simulations of the amplitude equa- 
tion JS]) and of the underlying equations of motion (|2]), with 
Tj = 0.02, h = 1, and all other parameters as in Fig. [1] Note 
that the simulations of the equations of motion, displayed on 
the right, show only the initial stage of the evolution that is 
simulated with the amplitude equation and displayed on the 
left (with T = 1 on the left equivalent to about m = 80 drive 
periods on the right, as discussed in the text), (a) and (b): 
k < kth — 1-36 and the soliton moves slightly and stops, (c) 
and (d): k > kth and the soliton splits into two. The so- 
lutions eventually settle into the known states of one or two 
stationary solitons. 



comes to a complete halt, as shown in Figs.[9fa) and (b). 

For certain parameter values we observe a noteworthy 
effect in which a boosted soliton splits into two. In order 
to estimate the threshold value kth for the wavenumber 
/c, above which this splitting occurs, we write the energy 
of the soliton using the Hamiltonian density that gives 
rise to the driven, but undamped, NLS equation. 



E 



ri 

J —CX) \ 



dtp 
dX 



+ I^Z-P - l^/'l^ + hRe{tP^) dX. (38) 



Following Eq. this energy evolves in time according 
to 



- |V^p + 2|V^|'^-/iRe(V^2) (39) 



It 



The right hand side of (|39|) is zero for = ae~*^sech(aX) 
and h cos(26>) = o? — 1 for any constant a — in particular, 
for the approximate solution V^app as well as the numeri- 
cally exact solution. By substituting the boosted approx- 
imate solution V^app(^)e"'^^^"^°^ into (|38|), we find its 
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FIG. 10: (Color online) A stable bound state of two solitons. 
The OS are absolute values of the displacements of the res- 
onators, obtained by a numerical integration of the discrete 
equations of motion Q , after a sufficiently long transient time 
has elapsed. The solid line is the stable solution obtained by 
solving the amplitude equation (|8]) as a boundary value prob- 
lem. The parameters are 7 = 0.565, /i = 1,7] = 0.3, and all 
others as in Fig. [T] 

energy to be 

/ 9x 2kiT{ai - 1) 2 o 
E{k) = 2a+ 1 + k^) + . + , - -al, 40 
smh(/c7r/a+) 3 ^ 

whereas for the static soliton E{k = 0) = 4a^/3. Thus, 
an obvious estimate for the threshold wavenumber kth re- 
quired to split a soliton into two is given by the condition 
E{kth) = 2E{0). For the parameters of Fig. El hh ^ 1-36, 
and indeed below this value the soliton does not split [in 
(a) and (b) k = 1.35] while above this value the soliton 
does split [in (c) and (d) k = 1.37]. At still larger values 
of k the soliton is destroyed by the boost and eventually 
decays to zero. For the parameters of Fig. [9] this hap- 
pens for k > 1.59. We note that although it might seem 
plausible to have values of k for which boosting a single 
soliton would split it into three, we were unable to detect 
such an effect. 



VIII. BOUND STATES 

We have considered the effects of pairwise interaction 
between solitons, and of boosting a single static soliton. 
It was shown by Barashenkov and Zemlyanaya [59] that a 
combination of both features within the framework of the 
PDNLS equation may lead to the formation of solitonic 
complexes, or bound states [g^. These complexes were 
found numerically, solving the PDNLS with an initial 
guess of the form 

V^b = V^(X-Xo)e^^(^-^°)+V^(X+Xo)e-^^(^+^°\ (41) 

where ipiX) = Asech{AX)e-'^ . For 7 = 0.565 and h = 
0.9 the rest of the parameters were found by means of a 



variational procedure elaborated in [59| to he = 6+, 
A = 1.14, k = -0.068, and Xq = 2.017. 

Using the PDNLS variational ansatz (|^T]) as an initial 
guess, we are able to obtain stationary solitonic bound 
states for the full amplitude equation ([8|) with 7^ > 0. 
Performing a linear stability analysis on these solutions, 
as described earlier using Eq. ([29|) , reveals that some of 
them are stable. Fig. [TOl shows one of these stable bound 
state solutions, obtained numerically using the amplitude 
equation ([8]), and nicely reproduced by a numerical inte- 
gration of the underlying equations of motion ([2]). 



IX. CONCLUSIONS 

We have investigated intrinsic localization of vibration 
in response to parametric excitation in an array of res- 
onators with a stiffening nonlinearity. Our analysis was 
chiefly performed on a single amplitude equation, which 
was derived directly from the underlying equations of mo- 
tion of the array to describe the slow spatio-temporal dy- 
namics of the system. The discreteness of the array im- 
poses an upper bound on the spectrum of linear modes. 
We have studied the case in which neighboring resonators 
oscillate out-of-phase, in the staggered mode, with an os- 
cillation frequency set slightly above the top frequency 
of the linear spectrum. One can similarly study ILMs 
in resonators with a softening nonlinearity, by changing 
the sign of in the equations of motion ([2]), and consid- 
ering the case in which neighboring resonators oscillate 
in-phase with an oscillation frequency set slightly below 
the bottom frequency of the linear spectrum. 

The array that we consider, hence also the amplitude 
equation we derive, are nonlinearly damped. Its localized 
modes emerge from two exact soliton solutions that ex- 
ist in the absence of nonlinear damping. We have shown 
that nonlinear damping increases the range of parame- 
ters for which localized solutions are stable. However, 
nonlinear damping increases the region in which the zero 
state is the only stable one, and it also stabilizes the non- 
zero uniform solution of the amplitude equation, which 
is modulationally- unstable if 7^ = 0. We have studied 
soliton interaction and soliton splitting, both in the pres- 
ence of nonlinear damping. We have also found a family 
of localized solutions in the form of bound states of two 
solitons. In a follow-up work, we intend to perform a 
more detailed investigation of the different localized so- 
lutions of the full amplitude equation (|8|) with 77 > 0, 
using numerical continuation. 

All results obtained from the amplitude equation are 
in excellent agreement with numerical solutions of the 
underlying equations of motion. This upholds the valid- 
ity of using a continuous PDE as a tool for analyzing 
ILMs, or discrete solitons, in a system whose original de- 
scription is given in terms of coupled ordinary differen- 
tial equations. Furthermore, our numerical simulations 
of the equations of motion suggest that the predicted ef- 
fects can be observed in parametrically-driven arrays of 
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real MEMS and NEMS resonators, thus motivating new 
experiments in these systems. 

We thank an anonymous referee for inquiring about 
the dynamical formation of solitons in our array, which 



prompted us to add Sec. |Vl This work was supported 
by the U.S. -Israel Binational Science Foundation (BSF) 
through Grant No. 2004339, and by the Israeli Ministry 
of Science and Technology. 
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